Melting of the classical bilayer Wigner crystal: influence of the lattice symmetry 
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The melting transition of the five different lattices of a 
bilayer crystal is studied using the Monte-Carlo technique. 
We found the surprising result that the square lattice has a 
substantial larger melting temperature as compared to the 
other lattice structures, which is a consequence of the specific 
topology of the temperature induced defects. A new melt- 
ing criterion is formulated which we show to be universal for 
bilayers as well as for single layer crystals. 

PACS numbers: 64.60.Cn, 64.70.Dv, 73.20.Dx 



Wigner crystallization of electrons on the surface of liq- 
uid helium was first observed experimentally by Grimes 
and Adams B. In the same year, Nelson, Halperin [Q, 
and Young developed a theory for a two stage con- 
tinuous melting of a two dimensional (2D) crystal which 
was based on the ideas of Berenzinskii [0] , Kosterlitz and 
Thouless H . Whether melting of a 2D crystal is a first or- 
der transition and proceeds discontinuously as described 
by the theories of Kleinert and Chui , or is a second 
order transition in which the crystal first transits into 
a hexatic phase retaining quasi-long-range orientational 
order and then melts into an isotropic fluid, is still an 
open question and a controversial issue. These studies of 
the melting transition of a 2D systems were directed to 
single layer crystals, which have the hexagonal symme- 
try. This is the most energetically favored structure for 
potentials of the form 1/r" ||. Disorder will influence 
Wigner crystallization as was demonstrated recently in 
Refs. @. 

In recent experiments on dusty plasmas JlO| | and on ion 
plasmas []TlJ| few layer and bilayer crystals were observed. 
Bilayer systems exhibit a much richer crystal structure 
(five different lattice types) as function of the inter-layer 
distance. This allows us to study the influence of the 
lattice symmetry on melting. Previously, the different 
types of lattices and structural transitions in a multilayer 
crystal at T = with parabolic confinement was analysed 
in Jl2|,|l3| . Different classes of lattices of the double-layer 
crystal were specified in |Q and in |15| the stability of 
the classical bilayer crystal was analysed in the harmonic 
approximation. 

In this letter we study the melting of a classical bilayer 
crystal, using the Monte Carlo (MC) simulation tech- 
nique. In the crystal phase the particles are arranged 
into two parallel layers in the (x, y)-plane which are a 
distance d apart in the z-direction. The layers contain 
equal density of particles n/2 and have close packed sym- 



metry. A single layer crystal is a limiting case of a bilayer 
crystal with d = and particle density n. 

We assume that the particles interact through an 
isotropic Coulomb (n — 0) or screened repulsive potential 
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where q is the particle charge, e the dielectric constant, 
r — (x, y,(z=0, d)) the position of the particle, and 
is the screening length. The type of lattice symmetry at 
T — depends on the dimcnsionless parameter v = d/ao, 
where d is the interlayer distance and clq = 1/ 
is the mean interparticle distance. For the classical 
Coulomb system (k = 0) there are two dimensionless pa- 
rameters v and r = q 2 /a^ksT which determine the state 
of the system. The classical Yukawa system (n ^ 0) at 
T =/= is characterised by three independent dimension- 
less parameters: v, T and A = na®. Below we measure 
the temperature in units of To = q 2 /a^ks and the energy 
in Eo = ksTo. 

The initial symmetry of the structure is set by the 
primitive vectors, the values of which are derived from a 
calculation of the minimal energy configuration for fixed 
v. In it was found that the bilayer Coulomb crys- 
tal exhibits five different types of lattices as function of 
the interlayer distance at T — 0: v < 0.006- hexago- 
nal, 0.006 < v < 0.262-rectangular, 0.262 < v < 0.621- 
square, 0.621 < v < 0.732-rhombic, and v > 0.732- 
hexagonal. Using the standard Metropolis algorithm |T^ ] 
we allow the system to approach its equilibrium state 
at some temperature T, after executing 10 4 -f- 5 x 10 5 
'MC steps'. Each MC step is formed by a random dis- 
placement of one particle. If the new configuration has 
a smaller energy it is accepted and if the new energy 
is larger the configuration is accepted with probability 
5 < exp(— AE/T), where S is a random number between 
and 1 and AE is the increment in the energy. In our 
numerical calculations the number of particles N may 
change for different types of bilayer crystals, but the par- 
ticle density remains the same. We took fragments of 
288 to 780 particles, where the shape of the specimen 
was determined by the T = crystal structure, and used 
periodic boundary conditions. Applying the Ewald tech- 
nique the potential energy is found by summation over 
all particles and their periodical images. 

The potential energy of the system as function of tem- 
perature is shown in Fig. 1. In the crystalline state the 
potential energy of the system rises linearly with tem- 
perature and then at some critical temperature it in- 
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creases very steeply. This denotes the beginning of melt- 
ing and is related to the unbinding of dislocation pairs, 
which we will discuss below. The square bilayer crys- 
tal [v = 0.4) exhibits a jump in the potential energy 
at melting of size 6 e = 0.71 x 10~ 2 /cbTo, and which 
is about a factor 2 larger than for a hexagonal lattice, 
i.e. at v = 0, S e = 0.39 x 10- 2 k B T , and at v = 0.8, 
S e = 0.31 x 10 _2 fcsT . Moreover, the square lattice 
has a substantial higher melting temperature, and conse- 
quently is more stable against thermal fluctuations than 
the hexagonal lattice. 

To characterize the order in the system we calculate 
the bond-angular order factor in each layer |17| 




and the translational order factor 

2 N/2 

Gl r = (^Y, e M*G-(n-r J ))), (3) 
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where index i — 1,2 refers to the top and the bottom lay- 
ers, respectively, and the total bond-angular order factor 
of the bilayer crystal is defined as Gg = (Gg + Gg)/2 and 
similar for G tr . N n b is the number of nearest neighbour 
particles (N n b = 6,4 for the hexagonal and square lat- 
tices, respectively), dj <n is the angle between some fixed 
axis and the vector which connects the jth particle and 
its nearest nth neighbour, and G is a reciprocal-lattice 
vector. 

From the behaviour of the order factors we can derive 
the temperature at which order is lost in the system. As 
seen from Fig. 2(a) the translational and orientational 
order is lost at the same temperature. Our numerical 
results show that for all five types of lattices the bond- 
angular order factor: 1) decreases linearly with increasing 
temperature (except very close to the melting temper- 
ature), and 2) it drops to zero just after it reaches the 
value 0.45. We found that Gg exhibits a universal behav- 
ior as shown in Fig. 2(b). We checked this for the bilayer 
crystal with screened and unscreened Coulomb interac- 
tion and for a single layer crystal with a Lennard-Jones 
V = 1/r 12 — 1/r 6 and a repulsive V = 1/r 12 interaction 
potential. From the present numerical results for Gg we 
formulate a new criterion for melting which we believe is 
universal: melting occurs when the bond-angle correlation 
factor becomes Gg w 0.45. 

Given the above mentioned criterion for melting we 
calculated the melting temperature using the harmonic 
approximation. Therefore, we numerically diagonalized 
the Hessian matrix for the finite fragment of the ideal 
structure at zero temperature of the crystal with period- 
ical boundary conditions in order to obtain the eigenval- 
ues. We checked that an increase of the size of the crystal 
fragment does not change our results. The melting tem- 
perature is then derived by linear extrapolating Gg to 



the value 0.45. In this way we obtained analytically T me i 
for different types of lattices which agrees with our MC 
calculations within 10%. 

Our results for the melting temperature are sum- 
marised in the phase diagram of Fig. 3 where we show 
the melting temperature as a function of v for two dif- 
ferent values of the screening parameter: A = for 
a Coulomb inter- particle interaction and A = 1 for a 
screened Coulomb interaction. For v = and A = we 
obtained the well-known value for the critical T = 132, 
resulting in T me i = 0.0076X"o- This critical value was first 
measured in Ref. [Q and found to be 137 ± 15. 

As seen in Fig. 3 the hexagonal (I and V) , rectangular 
(II) and rhombic (IV) lattices melt at almost the same 
temperatures. Further increasing the inter-layer distance 
we found that for v ~ 3 we obtained T me i w T me i(v = 
0)/v2- For the square bilayer crystal (phase III) the 
melting temperature increases up to T me / = 0.01078To 
with rising v and only for v > 0.4 we found that T me i 
starts to decrease with increasing v. It is surprising that 
the square lattice has a substantial larger melting tem- 
perature than the other lattices. This is true for Coulomb 
(A = 0) inter-particle interaction as well as for screened 
Coulomb. 

The detail analysis of the melting of the crystal in the 
vicinity of the structural phase boundary is much more 
complicated due to the softening of a phonon mode as 
shown in Ref. M| and is left for future work. To under- 
stand why the square lattice bilayer crystal has a consid- 
erable larger melting temperature, we investigated vari- 
ous temperature induced isomers of a single layer crystal 
and compared them with those of the square lattice bi- 
layer crystal with v = 0.4 which has the largest melting 
temperature. For bilayer crystals the topology of the de- 
fects is viewed as being composed by the top and the 
bottom staggered layers. Note, that the energy of the 
defects which occurs in the square lattice depends on 
the interlayer distance. At given temperature we found 
that during the MC simulation the system transits from 
one metastable state to another. They differ by the ap- 
pearance of isomers in the crystal structure which ap- 
pear with different probabilities. We found these iso- 
mers by freezing instant particle configurations during 
our MC steps. The topology of the defects, their energy 
and the bond-angular and the translational order fac- 
tors of these configurations are determined. Each point 
in Fig. 4(a,b) represents one configuration containing an 
isomer in a single layer and the square lattice bilayer 
crystals, respectively. The qualitative behaviour of both 
crystals during melting is similar although the energy 
of the defects in both lattice structures is substantially 
different. For the single layer (v — 0, Fig. 4(a)), all 
isomers at T\ = 0.00756T , just before melting, and at 
T 2 = 0.0076T just after melting, were obtained. Note, 
that for the square lattice (v = 0.4, Fig. 4(b)), we took 
Ti = 0.01076T and T 2 = 0.01078T . Typical calculated 
defect structures obtained from instant particle config- 
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urations freezed to T = arc shown in Fig. 5(a,b) for 
the hexagonal layer and in Fig. 5(c-f) for the square bi- 
layer crystal. First, at T = T\ the quartet of bound 
disclinations (Fig. 5(a)), point defects (Fig. 5(c,d)) and 
correlated dislocation (Fig. 5(e)) are formed. The point 
defects appear in pairs in our MC calculations, which are 
a consequence of the periodic boundary condition. Note 
that in a single layer crystal the total energy of a non 
bounded pair of a 'centred vacancy' and a 'centred in- 
terstitial' is E = 0.29fcsTo. In the square bilayer crystal 
the point defects like 'vacancy' and the 'interstitial', de- 
picted in Fig. 5(c,d), appear also in pairs and the energy 
of this unbounded pair is E = 0.315/cbTo. The discli- 
nations bound into a quartet and point defects produce 
only a negligible effect on the periodic lattice structure 
and G e = 0.8 + 0.9 and G tr = 0.85 0.95 (group A in 
Fig. 4(a,b)). It should be noted that in spite of prolonged 
annealing of the system during 5 x 10 5 MC steps at a tem- 
perature Ti, which is just below melting, we did not find 
more complex isomers than point defects and quartets of 
disclinations. 

At the temperature T = T 2 uncorrelated extended dis- 
locations with non-zero Burgers vector and unbounded 
disclination pairs are formed which causes a substantial 
decrease of the translational order (group B in Fig. 4(a,b) 
and defects shown in Fig. 5(b,f)). At this temperature 
single disclinations appear and the system looses order, 
both order factors become small and the system transits 
to the isotropic fluid (group C in Fig. 4(a,b)). 

Fig. 4(a,b) clearly illustrates that for a square lattice 
the defects which are able to destroy the translational and 
orientational order have a substantial larger energy than 
those of a single layer crystal with hexagonal symmetry. 
As a whole the localised and extended dislocations as 
well as disclinations in the square bilayer crystal are de- 
fects with a higher energy as compared to the ones in the 
hexagonal bilayer crystal. Thus, the square type bilayer 
crystal requires larger energies in order to create defects 
which are responsible for the loss of the bond-orientional 
and the translational order and thus for melting of the 
crystal. 

In conclusion, we studied the melting temperature of 
the five lattice structures in a bilayer crystal and found 
evidence that the melting temperature depends on the 
crystal symmetry. A square lattice has a substantial 
larger melting temperature than e.g. a hexagonal lat- 
tice. In order to understand this result we investigated 
the defect structures responsible for melting and found 
that the defects in a square lattice have a larger energy 
as compared to those in a hexagonal structure and conse- 
quently larger thermal energy is required to create them. 
We also formulated a new melting criterion: in two di- 
mensional layers and bilayers melting occurs when the 
bond-angular order factor is Go = 0.45, which is inde- 
pendent of the functional form of the interparticle inter- 
action. 
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FIGURES 

FIG. 1. The potential energy as function of temper- 
ature for the interlayer distances v = (solid circle), 
v = 0.4 (open squares); 

FIG. 2 (a) The bond-angular (Gg) and the transla- 
tional (Gtr) order factors as function of temperature for 
the interlayer distances v = (circles) and v = 0.4 
(squares); Gq: open symbols and G tr - solid symbols, (b) 
The bond angular order factor for different interaction 
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potentials: i) screened Coulomb: v = (solid squares- 
A = 1, open oncs-A = 3), and v = 0.4 (solid circles- 
A = 1, open ones-A = 3), ii) for the Lennard- Jones po- 
tential (solid rhombics), and iii) for the potential 1/r 12 
(open rhombics). 

FIG. 3. The phase diagram of the bilayer Coulomb 
crystal for without screening A — (open squares) and 
with screening A = 1 (circles). The vertical dotted lines 
delimit the different crystal structure which are depicted 
in the inserts (open symbols for the top layer and solid 
symbols for the bottom layer). The error bars denote 
the uncertainty in the temperature nearby the structural 
phase boundaries. 

FIG. 4. The bond-angular (solid squares) and the 
translational order factors (circles) of the different de- 
fects in (a) a single layer crystal [y = 0), and (b) in the 
square lattice bilayer system for v = 0.4. 

FIG. 5 The defects in a single layer crystal: (a) quar- 
tet of disclinations, and (b) two unbounded disclination 
pairs. In the square lattice bilayer crystal: (c) 'vacancy', 
(d) 'interstitial', (e) correlated dislocations and (f) a pair 
of disclinations. 
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